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Abstract. A characteristic particle method for the simulation of first order macroscopic traffic 
models on road networks is presented. The approach is based on the method particleclaw, which 
solves scalar one dimensional hyperbolic conservation laws exactly, except for a small error right 
around shocks. The method is generalized to nonlinear network flows, where particle approxi- 
mations on the edges are suitably coupled together at the network nodes. It is demonstrated in 
numerical examples that the resulting particle method can approximate traffic jams accurately, 
while only devoting a few degrees of freedom to each edge of the network. 



1. Introduction 

Macroscopic traffic models describe the evolution of the vehicular traffic density on a road by 
partial differential equations. A large network of highways, which includes ramps and intersections, 
can be modeled as a directed graph, whose edges are road segments that join and bifurcate at 
network nodes. The traffic density on each edge evolves according to some macroscopic traffic 
model, and at nodes, specific coupling conditions are specified. Since realistic road networks 
can easily consist of tens on thousands of edges, efficient numerical methods are crucial for the 
simulation, forecasting, now-casting, or optimization of traffic flow on networks. 

In this paper, we present a characteristic particle method for "first order" traffic models on 
road networks. The approach is based on the method particleclaw 130j . which solves scalar 
one dimensional hyperbolic conservation laws exactly, except for the immediate vicinity of shocks, 
where a small approximation error is observed. Particleclaw is used to evolve the numerical solution 
on each edge, and a methodology is presented for the coupling of the particle approximations on 
individual edges together through the network nodes. In recent years jTOJ, [T^l , particleclaw has 
been demonstrated to possess certain structural advantages over approaches that are based on a 
fixed grid. To name a few examples: 

• It possesses no numerical viscosity, and thus preserves a reasonable accuracy even with 
very few degrees of freedom. This is in contrast to Godunov's method [15] and other low 
order schemes. 

• It is optimally local, in the sense that particles move independently of each other, and 
communication between neighboring particles occurs only at shocks. In contrast, high 
order finite difference approaches [22], finite volume methods [31] . or ENO [T7J/WENO 
[26] schemes use wide stencils to achieve high order, which poses a particular challenge 
right at network nodes. 

• It is total variation diminishing, yet it is second order accurate even in the presence of 
shocks [TT]. In contrast, many limiters in fixed grid approaches need to reduce the order 
of convergence near shocks in order to avoid overshoots. 
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• The approximate solution is represented as a piecewise similarity solution of the underlying 
equation, that is continuous except for right at actual shocks. In contrast, the reconstruc- 
tion in finite volume [3T] and Discontinuous Galerkin [2E1 E] methods tends to possess 
spurious discontinuities at every boundary between cells. 

• It is adaptive by construction. Particles adapt to the shape of the solution, and the 
approach can be generalized to incorporate stiff source terms [T3]- This is in contrast to 
fixed grid methods, which would need to use mesh refinement techniques [3. for adaptivity. 

These advantages render particleclaw as an interesting candidate for the simulation of traffic 
flow on networks. The need for the design of specialized numerical schemes for traffic models 
on networks has been pointed out for instance by Brctti, Natalini, and Piccoli j4], who design 
an efficient numerical approach for a very specific class of flux functions, for which it suffices to 
track the location of the transition between free and congested traffic flow. The idea of tracking 
features is also shared by front tracking methods [21], which approximate the whole solution by 
finitely many moving discontinuities. In fact, in analogy to front tracking, particleclaw can also 
be interpreted as a rarefaction tracking |12j . 

This paper is organized as follows. In Sect. [2j the class of traffic models under consideration 
is outlined, and the coupling conditions for road networks are given. The characteristic particle 
method particleclaw is presented in Sect. [3j and in Sect. [4j we demonstrate how the approach 
can be generalized to nonlinear flows on networks. Numerical results are shown in Sect. [5] and in 
Sect.[6j we present conclusions and give an outlook. 

2. Traffic Models on Highway Networks 

2.1. Macroscopic Traffic Models. Macroscopic traffic models treat the traffic as a continuum, 
and use partial differential equations to describe the temporal evolution of the vehicle density 
u(x, where x is the position on the road (the flow is averaged over all lanes that go in one 
direction), and t is time. If vehicles move with the velocity v(x, t), then the conservation of vehicles 
(in the absence of ramps) is described by the continuity equation u t + (uv) x = 0. The assumption of 
a direct functional relationship between the velocity and the density, v = v(u), yields the classical 
Lighthill-Whitham-Richards (LWR) model [23 [2PJ 

«t + (/(«))* = , (1) 
which is a scalar hyperbolic conservation law with a flux function f = uv that equals the traffic 
flow rate. Due to the direct density-velocity relationship, the LWR model does not model vehicle 
accelerations, and thus is not able to describe non-equilibrium features such as phantom traffic 
jams or self-sustained traffic waves ("jamitons") |14j . To overcome this limitation, numerous types 
of "second order" models have been developed, such as the Payne- Whitham model [27], the Zhang- 
Aw-Rascle model (32[[2], phase transition models [8], and many more. These introduce a velocity 
(or velocity-like) variable as an independent unknown into the equations, resulting in a system of 
balance laws. In this paper, we restrict our study to effects that can be captured reasonably well 
by the scalar LWR model, i.e. we arc interested in the large scale, nonlinear equilibrium behavior 
of traffic flow. 

As it is common in traffic models [TS], we assume that the LWR flux function satisfies the 
following conditions: (i) no flow for empty road and bumper-to-bumper densities, i.e. /(0) = 
= f(u m ), for some maximum density u m ; and (ii) concavity, i.e. f"(u) < Vu € [0,it m ]. As a 
consequence of these assumptions, there is a unique maximum flow rate /* = f(u*) that occurs for 
an optimal density u* . Two parameters of / are fundamentally important: the slope v m = /'(0) 
is the velocity of vehicles when alone on the road (i.e. approximately the speed limit); and u m 



Even though densities are commonly denoted by p, here we use u, in order to express the fact that the numerical 
approach applies to more general network flows. 
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is the number of lanes divided by the average vehicle length plus a safety distance. Other than 
for these two values, the precise functional shape of / can depend on many factors, such as the 
road geometry, the drivers' state of mind, etc. Frequently, one assumes a simple parabolic profile 
f(u) = v m u(l — This particular shape was inspired by measurements of Greenshields |16j . 

and even though it does not fit well with contemporary measurements |18] , it is commonly used 
due to its simplicity. We shall do so here as well, as it simplifies the presentation of the numerical 
approach in Sect. [4j 

2.2. Traffic Networks. A traffic network is a directed graph of network edges (road segments) 
that join and bifurcate at network nodes. On each edge, the traffic evolves according to the LWR 
model ([!]) with a flux function that is specific to this edge (here, we assume that on each edge the 
flux function does not explicitly depend on space or time). In order to have a feasible model for 
traffic flow on road networks, one must formulate coupling conditions that guarantee the existence 
of a unique entropy solution, given suitable initial and boundary conditions. This problem was 
first addressed by Holden and Risebro [22] , and subsequently generalized by various other authors 
(see below). In the following, we briefly outline the key ideas of coupling conditions at network 
nodes. 

At a network node, one has the following setup. Let the node have n ingoing edges (roads), 
numbered i = 1,2, . . . , n, each of which carries a vehicle density Ui(xi, t), where Xi G [0, Li], and a 
flux function /j (u) . Similarly, let the node have m outgoing edges, numbered j — n + 1, . . . , n+ m, 
each carrying a density Uj(xj,t), where Xj G [0, Lj], and a flux function fj(u). Ingoing edges end 
at the node at x% — Li, and outgoing edges start at the node at Xj — 0. The conservation of 
vehicles requires that the total inflow into the node equals the total outflow out of the node, i.e. 

n n-\-m 

£/ i (« i (L i ,t))= Yl /i(%(<M))Vt. (2) 

i—l j=n-\-X 

In order to obtain further rules for the temporal evolution of the solution at nodes, a generalized 
Riemann problem is formulated, as follows. Let all edges be extended away from the node to ±oo, 
and consider initial conditions that are constant on each edge, i.e. 0) = it,-, for x G [-co, Lj], 
and Uj(x,0) = Uj for x G [0, oo]. The question is then: what are the new states Ui — Ui(Li,t) 
and Uj = Uj(0,t) at the node for t > 0? These new states define the solution at all times, since 
the problem admits a self-similar solution Ui(x,t) = Wi( 7 X ) and uj(x,t) = Wj(j). Let the flux 
values of the new states be denoted 7$ = fi(ui) Vi G {1, . . . ,n + m}- Clearly, these fluxes must 
satisfy the conservation condition ([2]), i.e. 

n n+m 
i—l j—n+1 

Moreover, the new states Ui Vi € {1, . . . , n + m} must generate a solution for which information 
on each edge is either stationary or moves away from the node (information cannot move "into" 
the node). As shown in [22], this information direction requirement implies that the new fluxes 
satisfy the inequality constraints 

7i G SI; , where fij = [0, /j(min{tt,, u*})] Vi G {1, ... ,71} , and (4) 
7j G ftj , where Qj = [0, fj(max{uj, u*})] V j G {n + 1, . . . , n + m} . (5) 

In words: on ingoing edges, the new flux must be less than or equal to the demand flux, which 
is given by the old flux if Ui < u*, and the maximum flux /* if m > u*. Likewise, on outgoing 
edges, the new flux must be less than or equal to the supply flux, which is given by the old flux if 
u j _^ u* , and the maximum flux f* if Uj < u* . 

Conditions Q, and |5]) allow infinitely many possible new fluxes. Further conditions must 
be provided to single out a unique solution. A reasonable set of additional conditions is given by the 
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drivers' desired destinations matrix, denned as follows. Let 7i n denote the vector of ingoing fluxes, 
and 7 ou t be the vector of outgoing fluxes. The fact that very specific percentages of drivers that 
enter the node through edge i exit the node on edge j (different drivers have different destinations) 
can be encoded via the linear system 

A ■ 7i„ = 7out , (6) 



where 



0-l,n+l • . . 0-n,n+l 



hi 



A 



Tout = 



The matrix A is column-stochastic, i.e. all dij € [0, 1] and X^="+i a i,j ~ Thus, the flux balance 
condition ^ is automatically guaranteed by relation ([6|, since ■ 7* ou t = • A ■ 7; n = ■ 7i n , 
where e = (I, . . . , I) T . Condition ^ together with the constraints Q and ^ yields that 7i n must 
lie in the feasibility domain 

= {7 e Ox x ■ • • x Q n I A ■ 7 e O n+1 x • • • x f2„+ m } c R™ . (7) 

The selection of a specific 7i n G Q requires further modeling arguments. Possible criteria are, for 
instance, entropy arguments [52], the modeling of the intersection geometry [T!5], or simply the 
assumption that drivers behave such that the throughput through the node is maximized [5] , 
i.e. one solves the linear program 

max e T ■ 7i n s.t. 7; n € Q . (8) 

In the examples presented in Sect. |5j we shall follow the latter option, even though the other 
alternatives are possible as well. The modeling has to be augmented by one small but important 
detail. It is possible that (JsJ) does not possess a unique solution, namely if the extremal boundary 
of f2 is perpendicular to e. In this case, one can introduce the additional constraint that 7i n € cR, 
where c £ WL n is a given constant that models the merging behavior at the node. 

The definition of the generalized Riemann problem is finalized by selecting new states as follows. 
On ingoing roads, choose Hi = Ui if 7; = fi{ui), otherwise choose u, > u*, s.t. fi(v.i) — 7,. 
Similarly, on outgoing roads, choose Uj = Uj if jj — fj(uj), otherwise choose Uj < u*, s.t. fj(uj) = 
7j. By construction of ^ and any resulting shocks and rarefaction waves are guaranteed to 
move away from the node (i.e. forward on outgoing and backward in ingoing edges). 



3. Particleclaw 



3.1. Characteristic Particles and Similarity Interpolant. On each network edge, we have 
a scalar one-dimensional hyperbolic conservation law 

«t + (/(«))* = , u(a: > 0)=«o(a:) ) (9) 
where the flux function / is assumed to be twice differentiable and concave (/" < 0) on the 
range of function values (see [TT] for extensions of the approach to flux functions with inflection 
points). We consider a special subset of exact solutions of ([9|, which can be represented by a finite 
number of characteristic particles, as follows. A particle is a computational node that carries a 
position x l (t), and a function value u l (t). Note that we shall consistently denote particle indices 
by superscripts, while subscripts are reserved for edge indices. In the following, for convenience, 
we omit the time-dependence in the notation. Consider a time-dependent set of n particles P = 
{(x 1 , it 1 ), . . . , (x n , «")}, where x 1 < ■■■ < x n . On the interval [a; 1 , a;"] that is spanned by the 
particles, we define the similarity interpolant Up{x) piecewise between neighboring particles as a 
true similarity solution of (J9J) . If u l 7^ u l+1 , it is implicitly given (and uniquely defined, since / is 
concave) by 

s-s* ^ f(Up(x))-f'(u l ) fm 
x i+i _ x i f'(u i + 1 ) - /'(u*) ' 1 ' 
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If u % = u t+1 , the interpolant is simply constant Up{x) = u % . As shown in [TT], the interpolant 
Up is an analytical solution of the conservation law ^ , as each particle moves according to the 
characteristic equations 

x = /» , u = , (11) 
i.e. {x l {t),u l {t)) = {x l {0) + /'(u l (0))t,?i t (0)). The reason for this fact is that due to the par- 
ticular form of Up, each point (x(t),u(t)) on it does move according to the same characteristic 



equations (11). Note that strictly speaking, Up(x,t) is a solution only in the weak sense, since 



the derivative ^Up(x,t) is discontinuous at the particles. However, since Up is continuous, the 
Rankine-Hugoniot shock conditions [5] are trivially satisfied. 

The interpolant Up is called "similarity interpolant", since the solution between neighboring 



particles (10) is either a rarefaction wave that comes from a discontinuity (if f'(x l ) < f'(x l+1 )) 
or a compression wave that will become a shock (if f'(x l ) > f'(x l+1 )). As a consequence, the 
approach can be interpreted as rarefaction tracking |12j . This expresses both its similarities and 
differences to front tracking approaches |21l 123] , which approximate the true solution by a finite 
number of shocks. 

Just as the true solution of Q may cease to be continuous after some critical time, the similarity 
approximation exists only up to the time of the first collision of neighboring particles. For a pair 
of neighboring particles (x l ,u l ) and (x l+1 ,u t+1 ), the time of collision (i.e. they have the same 
cc-position) is given by 

T l = - - (12) 

-/'(«*) ' 

given that f'(x l ) > f'(x l+1 ). Consequently, for n particles, the time of the first collision is 
T* = min ({T l | T % > 0} U oo). At that time, a shock occurs (at x i = x i+1 , from u i+1 to u*), and 
the method of characteristics cannot be used to evolve the solution further in time. 



3.2. Representation of Shocks. The characteristic motion of particles, as described in Sect. |3.1[ 
can only be performed if no shocks are present in the numerical solution. One idea to overcome 
this limitation, thus admitting an evolution of solutions with shocks, is to merge particles upon 
collision. This approach was first presented in [10] . analyzed and generalized in [TH [12], and 
implemented in the software particleclaw [30] . The merging of two particles (x l , it*) and {x l+1 , u i+1 ) 
with x % — x t+1 into a single new particle (x l ,u l ) is performed such that the total area under the 
similarity interpolant Up is exactly preserved. As shown in |llj . the area under Up between two 
neighboring particles equals 

U P (x) dx = {x l+1 - x 1 ) a{u\u l+1 ) , (13) 
where 

a(ttW » )= ir<Mhz01 (14) 

[f'W] ai 

is a nonlinear average function (note that [f(u)]™ — f(w) — /(«)). Equating the area under the 
interpolant before and after the merge, in general, yields a nonlinear equation for u l , which can 
be solved up to machine accuracy by a few Newton iteration steps (the geometry of the problem 
generally yields a very good initial guess) . In the case we are solving here, the interpolant is linear 
which simplifies many calculations. 

The merging of colliding particles replaces a discontinuity by a continuous interpolant, and thus 
the numerical approximation can be evolved further in time using the characteristic particle motion 
(11). Since the merging of particles i and i + 1 modifies the interpolant in the interval [x 1- 1 , x l+2 ] , 



this approach introduces a small error right around shocks. We control the magnitude of this error 
by the following additional step. Let a parameter d be given on the edge that provides an upper 
bound on the distance of particles adjacent to a pair of particles that need to be merged. Consider 
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two particles i and i + 1 that need to be merged, because x l = x l+1 and f'{u l ) > f'(u l+1 ). 
If x % — x 1 " 1 > d, then a new particle is inserted (on the interpolant) at x l — d. Moreover, if 
x x+2 — x x+1 > d, then a new particle is inserted at x l+1 + d. After the relabeling of the particle 
indices, the original particles i and i + 1 are merged, in the way described above. 

Note that there is an alternative to particle merges: the generalization to shock particles, 
as introduced in [T3]. The advantage of this version of the particle method is that shocks are 



represented exactly. The price to pay is a particle evolution that is more complicated than (111. 
In this paper, for the application to macroscopic traffic models on networks, we consider the 
approach that uses characteristic particles. While it is certainly possible to formulate the ideas 
presented below with shock particles, the simplicity of the characteristic particle motion admits 
an easier presentation and analysis of the methodology. Moreover, the approach presented here 
introduces an intrinsic error at network vertices, and thus the exact nature of shock particles is 
less of an advantage. 

4. Generalizing Particleclaw to Network Flows 

In this section we demonstrate how the approach particleclaw can be generalized to nonlinear 
traffic models on highway networks. Our goal is to obtain an overall approach that is highly 
modular, in the sense that the solution on each edge can be evolved independently of the other 
edges, and edges communicate only during a synchronization step. The specific methodology is as 
follows. Having sampled the initial conditions on each edge and synchronized them (as described 
below), we pick a time step At. During each time increment t £ [t n ,t n + At], the solution on each 
edge is evolved using the simple particle method described in Sect. [3j with a special treatment at 



the first and the last particle (see Sect. 4.1 ). At the end of the time step, at each network node, the 
adjacent edges are synchronized with each other: first, the numerical solution (which may have 
partially moved away from the edge) is interpolated/extrapolated back onto the edge; second, 



the generalized Riemann problems (described in Sect. 2.2) are invoked; lastly area is suitably 
re-distributed. All these operations are performed such that the approach is exactly conservative, 
i.e. no area is lost or created. 

The independent evolution of the solution on the edges renders the approach extremely adapt 
towards parallelization: each edge can be stored in its own share of memory, communication 
between the different edges need only occur during the synchronization, and very little information 



must be transferred (see Sect. 4.2). This methodology is possible due to the finite speed of 
information propagation of the hyperbolic conservation law ([I]) on each edge. There is a maximum 
synchronization time step Ai max , such that for all At < Ai max , information does not propagate 
further than half the length of each edge. As a consequence, in the synchronization step, all 
nodes can be treated independently of each other. Note that the maximum admissible time step 
At m&x between synchronization events is on the order of the smallest edge length divided by the 
fastest characteristic velocity. This is significantly larger than the maximum admissible time step 
in many traditional numerical approaches, which is on the order of the grid resolution divided 
by the fastest characteristic velocity. Hence with the presented particle approach, the relatively 
costly generalized Riemann problems need to be called much less frequently. 

It should be pointed out that even though the method is exactly conservative, it is not exact. 
Due to the uncoupled evolution of the edges, a certain amount of information is lost, resulting in 
approximation errors that increase with the size of At. Thus, there could be accuracy constraints 
on At that are more stringent than the stability constraints. Below, we outline the required 



additions to particleclaw (Sect. |4~T[), describe the synchronization step (Sect. 4.2), and show how 



exact area balance is achieved (Sect. 4.3) 



4.1. Virtual Domain, Excess Area, Virtual Area, and Area Credit. One key idea of the 
particle approach is that the numerical solution can be advanced on each edge, without considering 
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the coupling with other edges at the network nodes. Clearly, this evolution incurs an error near 

we 



the two edge boundaries, since the coupling information is not used. However, in Sect. 4.2 
design the coupling step such that the arising gaps in information near the edge boundaries can 
be filled (with satisfactory accuracy) during the coupling step. 

The coupling-free evolution on an edge x G [0, L] works very similarly to the basic particleclaw 
methodology described in Sect.[3j with a special twist at the first and the last particle, as described 
below. When no particles share the same x-position, then all particles are simply moved according 
to the method of characteristics pdj ). Since this applies also to the first and the last particle, 
the domain of representation of the solution [x 1 (t) , x N (t)] does in general not match with the 
computational domain [0, L]. In the case that particles extend beyond the edge into the virtual 



domain, i.e. x < or x > L, the similarity interpolant ( 10 ) defines the numerical solution in 



particular up to x = or x = L, respectively. In addition, it defines an excess area j 1 Up(x) dx or 

x N 

J L Up{x) dec, respectively. In the case that particles do not reach the edge boundary, i.e. x 1 > 
or x N < L, we can define an extrapolation on [OjX 1 ] or [x , L], respectively, if we are given a 

value for the virtual area u Ax or f^ N u dx, respectively. How to construct this extrapolation 
is described in Sect. 14.31 



The merging of colliding particles works as described in Sect. |3.2[ with one exception. If for the 
merging of particles i and i + there is no particle i — then a new particle is added at x % — d with 
the same value as u l , and a "left-sided" area credit II = du l is recorded. Similarly, if there is no 
particle i + 2, a new particle is added at x i+1 + d with the same value as and a "right-sided" 

area credit In = du l+1 is recorded. After these insertions, the merge can be performed. 

4.2. Synchronization Step. The synchronization of area between the edges happens by moving 
area from edges in which information is flowing into the node to edges from which information is 
flowing out of the node. As introduced in Sect. |2.2[ we denote the edges on which vehicles enter 
the node "ingoing". On the ingoing edge i G {1, . . . , n}, the position of the node is at x — Li. 
Conversely, the edges on which vehicles exit the node are called "outgoing" . On the outgoing edge 
j G {n + 1, . . . , n + m}, the position of the node is at x = 0. Since the LWR model ([T]) admits 
characteristic velocities of either sign, information can propagate either into or away from the 
node, both for ingoing and outgoing edges. Therefore, we introduce further terminology: edges for 
which information is going into the node are called "influencing" , edges for which information is 
going away from the node are denoted "affected" , and edges with zero characteristic velocity are 
called "neutral" . When solving a generalized Riemann problem at a node, the procedure described 



in the last paragraph of Sect. 2.2 implies that all edges for which the flux changes (from fi(ut) to 
7i) become automatically affected or neutral. Influencing edges arise if % = fiiui) and in addition: 
itj < u* for ingoing edges and Ui > u* for outgoing edges. 

Let us now focus on one node in the network, with n ingoing and m outgoing edges. In the 
following, subscripts denote the edge index, and superscripts denote the particle index. With this 
notation, the last particle on the ingoing edge i G {1, . . . ,n} is (x^',uf z ), and the first particle 
on the outgoing edge j G {n + 1, . . . ,n + m} is (xj, uj). We design our approach such that at the 
end of the synchronization step (and thus at the beginning of the evolution on each edge) the first 
and last particle on each edge arise as solutions of the generalized Riemann problems described in 



Sect. 2.2 Consequently, due to ([6]), we have the m conditions 



\ " t , N 



/_^ a i,jfi{i J 'i *) = fj( u j) Vj€{n + l,...,n + m}, (15) 



given by the desired destinations matrix A, and — if necessary — further n — 1 conditions 

A(uf') = V»€{l,...,n} (16) 
for some j3 G K, given by the merging vector c G R". 
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We no w de rive, first for ingoing edges, an evolution equation for the excess area or virtual area 



(see Sect. 4.1) that is between the end of the edge Li and the last particle x i i . Letting u denote 



a true solution of the hyperbolic conservation law Q, this area Ii = f£* u(x,t) dx has the rate 
of change 



(17) 



j t h = fi(ui{L h t)) - /i(uf') + uf*/i(wf 4 ) ■ 

Here, the first two terms come from the PDE ([I]) , and the third term stems from the fact that . 
moves with velocity fKu?*). Analogously, we compute the evolution of the excess/virtual area 
Ij = f^ 1 u(x, t) dx on an outgoing edge as 



dt Ij 



/i(%(0.t))-/j(«J) + «j^(«i) 



(18) 



At the beginning of the coupling-free evolution, the first and last particle are on the edge 
boundaries, and thus all excess/virtual areas are zero. At the end of the coupling-free evolution, 
at the considered node, we know the excess area for all influencing edges. Moreover, for all neutral 
edges the excess area is zero. The idea is now to use the evolution equations ( 17 ) and ( 18 1 , together 



with the relations ( 15 ) and ( 16 ), to determine the unknown virtual areas at all affected edges. Any 



area credit that has been taken due to merges at the first or last particle one an edge (see Sect. 4.1 ) 
is simply subtracted from this balance. Below, we present this idea for three fundamental cases: 
a bottleneck, a bifurcation, and a confluence. Many practically relevant road networks can be 
constructed from these three types of nodes. 



4.2.1. Bottleneck. A bottleneck is a sudden change of the road conditions, i.e. n = 1 and m = 1. 
Here, the feasibility domain Q is fl = [0,dx] n [0, S2], where d\ = /i(min{iii, u\}) and s 2 = 
/2(max{w2, Thus, the solution of the optimization problem ^ is 71 = min{<ii,S2}, which 

allows for two cases. If 71 = d±, then all vehicles can pass through the node, and consequently 
edge 1 is influencing (or neutral), and edge 2 is affected (or neutral). Otherwise, if 71 < d\, then 
not all vehicles can pass, and a jam is triggered. Consequently, edge 1 is affected, and edge 2 is 
influencing or neutral. 

In either case, since the influx equals the outflux at all times, we have that fi{u±{L±, t)) — 
/2 (1*2(0, t)) and fiiu^ 1 ) = j '2(^2), and thus the evolution equations (17) and (18) imply that 

d 



dt 



which is a known time-independent quantity. Together with I\(t) = I 2 (t) = 0, we obtain a relation 
I 2 {t + At) — Ii(t + At) = CAt that allows to obtain I 2 from I\ (in the case 71 = d\), or 1\ from 
I 2 (in the case 71 < d\). Figure [l] shows an example of the movement of particles into the virtual 
domain, and the evolution of excess area (Jx) and virtual area (^2)- 



4.2.2. Bifurcation. In a bifurcation, one road splits into two, i.e. n — 1 and m — 2. The relative 
numbers of vehicles that exit the node on edge 2 or edge 3 are given by ai j2 and 01,3 = 1 — 
ai,2, respectively. Here, the feasibility domain |7| is Q = {71 e [0, di] \ aijji € [0, Sj] Vj € 
{2,3}, where d\ = /i(min{iti, u\}) and Sj — fj(max{uj,u*}) Vj G {2,3}. The solution of the 
optimization problem Q is given by 71 = min{<ii, ^^'^jl- This allows two possibilities. If 
71 = dx, then all vehicles can pass through the node, and edge 1 is influencing (or neutral), and 
edges 2 and 3 are affected (or neutral). Otherwise, if 71 < d\, then not all vehicles can pass, and 
edge 1 becomes affected. Of the outgoing edges, the one that causes the congestion is influencing 
or neutral, while the other outgoing edge is affected or neutral. 



In either situation, we can use the evolution equations (17) and (18) to determine the virtual 
areas, since ( 15 ) provides two conditions that relate I\, I 2 , and I3 with each other. As an example, 
consider the case of edge 2 being influencing, and edges 1 and 3 affected. Here, I 2 is known (by 
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Figure 1. Area tracking for a bottleneck. The excess area 7i allows the computation of the virtual area 
Ii. 



integrating the interpolant ( 10 ) ) . From (17) and (18 1 we first compute C = — = 

"f'/iOf 1 ) - ^ulfeiul), and then obtain h(t + At) = ^h(t + At) + CAt. Knowing I u we 
can then determine I3, by considering ^(^3 — ai,3h)- 

4.2.3. Confluence. Here, two roads converge into one, i.e. n = 2 and m = 1. Flux conservation 
simply states that 71+72 = 73- The feasibility domain |7]) is now two-dimensional: ft = {(71, 72) € 
[0,di] x [0,d 2 ] I 7i+72 € [0, s 3 ]}, where di = ^(min-t^, u*}) \/i e {1, 2} and s 3 = / 3 (max{u 3 ,^}). 
For the optimization problem (JsJ) , we distinguish two cases. If d\ + c?2 < S3, then all vehicles can 
pass through the node, and the ingoing edges are influencing (or neutral), and the outgoing edge 
is affected (or neutral). In contrast, if d\ + c?2 > S3, then a jam is triggered. How the backwards 
going information distributes among the two ingoing edges depends now on the merging vector 
(01,02). The following outcomes are possible: (i) all vehicles from edge 1 get through, i.e. 71 = di 
and a jam occurs on edge 2, i.e. 72 < di. In this case, edge 1 is influencing, edge 2 is affected, and 
edge 3 is influencing or neutral; (ii) a jam occurs on both ingoing edges, i.e. 71 < d\ and 72 < g?2- 
Here, both ingoing edges are affected, and the outgoing edge is influencing or neutral; (iii) is the 
same as (i), but with edges 1 and 2 reversed. 

Again, the virtual areas at the affected edges can be obtained systematically from the known 
virtual areas at the influencing and neutral edges. In the case that all vehicles pass through the 

node, we have I 3 (t + At) =h(t+ At) + J 3 (t + At) + (ulf^ul) - < 1 f[ (uf 1 ) - u^V^^ 2 )) A*- 
In turn, if jamming occurs, e.g. on both ingoing edges, we use the knowledge of the merging ratios 
3^ = ^ to determine I\ and I2 from I3. The methodology is analogous to the cases presented 
above. 



4.3. Representation of Virtual Area by Particles. In Sect. |4.2( we have described how the 
virtual area at affected edges can be determined. The presented methodology tracks area exactly, 
and thus the numerical approach is exactly conservative. In order to finalize the synchronization 
step, we must account for the imbalance of area by changing parts of the solution, thereby creating 
or removing area. 

Throughout the design of particleclaw, the two main guiding principles are: (i) conserve J u dx 
exactly and as locally as possible; and (ii) ensure that the approach is TVD, i.e. no spurious 
overshoots or undershoots are generated. We conclude the synchronization step by modifying the 
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numerical solution, while following these two principles. The conservation principle means that 
the excess area and the virtual area must be converted into actual area under particles on the edges 
near where that area was recorded. The locality requirement is respected by a finite domain of 
influence: the solution is never modified further than At times the maximal particle velocity from 
the node. The TVD condition means that the modified solution does not possess values above or 
below the range of values that is presented in the nearby particles and the solution u that comes 
from the generalized Riemann problem ^ . The implementation outlined below violates the TVD 
condition only under exceptional circumstances. In fact, in the numerical experiments presented 
in Sect.[5j this non-TVD "last resort" solution turned out to never be required. 

The precise methodology is as follows. First, the numerical solution is truncated or extrapolated 
to the edges, by adding a particle either on the interpolant, or with the same value as the closest 
particle, respectively. All excess areas (i.e. the parts that extend beyond the domain of each edge) 
can now be removed, since their contribution has been accounted for in the step presented in 
Sect. |4.2[ Now, the generalized Riemann solver is called with the new values on the node position. 
The resulting w-values are added to the solution as new extremal particles (unless the value has 
not changed, u — u, in which case no particle is added) with the same x-value as the extremal 
particle. Being able to have two particles with the same x-value is a peculiarity of the particle 
method, and it implies that the particle insertion does not change the area under the solution. 

With this, on each edge we have a numerical solution, whose extremal states are solutions of 
the generalized Riemann problems. What is left is the accounting for the correct balance of area. 
We do so by modifying the numerical solution, however, without changing the previously found 
new extremal states. This guarantees that at the end of the synchronization step, the fluxes at 
the all the nodes are solutions of the generalized Riemann problem. In the preceding subsections 
we have derived the area I that is supposed to be under the solution between the network node 
and the extremal particle. In general, the area under the numerical solution constructed above by 
simple interpolation/extrapolation has a small disparity with /. The difference, A/, needs to be 
added to the affected edge, by modifying the numerical solution. Without loss of generality let 
us examine an outgoing affected edge. We look for a new solution to replace the current solution 
near the end of the edge. This new solution must have an area difference AI from the area of the 
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current solution, and the first particle must remain unchanged in order to satisfy the generalized 
Riemann problem. We follow the steps below (starting with k = 2) to represent the required 
change in area using particles (illustrated in Fig.[2|: 

(1) Attempt a solution that comprises a constant part and a linear part so that it connects 
particles 1 and k in a continuous fashion. If such a solution can provide AI, accept it. 

(2) Attempt a solution with constant value u in the domain [xi,Xfc], where u € [mini<fc(itj), 
maxj<fc(tii)]. If such a solution can provide AI, accept it. If this solution is not accepted, 
and Xk < At max u (/'(u)), increase k by one and go back to (1), otherwise, go to (3). 

(3) Use a solution that has a constant value u in the domain \x\ , x^] with the value of u needed 
to match AI. 

In the algorithm above, step (3) is a fail-safe mechanism that will always result in a solution, but 
the solution may be non-TVD or even unphysical. 

Figure [2] shows eight exemplary cases using the same initial particle configuration but different 
AI. In these examples, k could not increase beyond 4: AI = 0.16 resulted in a solution at step 
(3) with k = 4; AI = 0.08, 0.04, 0.021, with step (1) and k = 4; AI = 0.015, 0.005 with step (1) 
and k = 3; AI = -0.025 with step (2) and k = 4; and AI = -0.065 with step (3) and k = 4. 

Performing this area reconstruction on all affected edges guarantees that area is preserved 
exactly; it may increase or decrease the number of particles. 

5. Numerical Results 

5.1. Bottleneck Test Case. As a first test case, we consider the simplest non-trivial network: 
a bottleneck that consists of two edges, each of length L\ = L2 = 1 that are joined linearly, 
as shown in Fig. [3j In dimensionless variables, the maximum traffic densities on the two roads 
are it™ = 2 and u'f = 1, and the maximum vehicle velocities are = 1 and v 2 n — 1-5- This 
example can be interpreted as a model for the situation of two lanes merging into one lane, and 
at the same position the speed limit increasing by 50%. The initial conditions are U\{x) = 1 — x 
and 1*2(2;) = 0.8 x, and the boundary conditions of Ui(0, t) — 1 and 11,2(1, t) = 0.8 are prescribed 
whenever information is entering the edge through these boundaries. The final time is ta na \ = 3. 
In Fig. [3] (and also in Fig. [5]), the solution is shown in two ways: the shade of gray becomes darker 
with higher values of u, and the plot of u vs. x is overlaid with dots representing the particles of 
the method. Each edge is annotated with the maximum velocity v m and two arrows showing the 
direction of vehicles. The thickness of each segment is proportional to its maximum density u m . 

It is shown in |llj that particleclaw itself is second-order accurate with respect to the maximum 
distance of particles (away from shocks or when using shock sharpening). What we investigate 
here is the error as a function of the size of the synchronization time step At, while having a large 
number of particles, such that the error due to the spacial approximation is negligible. Specifically, 
we provide the following two parameters for particleclaw: an initial/desired distance of particles 
of h = 8 • 10~ 5 is given, and the distance of particles that are inserted near a shock is d = 2 • 10~ 5 . 
Both parameters result in a very fine resolution, using more than ten thousand particles on each 
edge. We choose the time step from At — 2~ k / 2 , where k € {4, . . . , 15}, and compute a reference 
solution using k = 20. Since here we do not wish to measure the error at a shock, we evaluate 
the error in the segment x € [0,0.3] on edge 2, in which the solution is smooth, but non-linear. 
Figure [4] shows the L°°([0, 0.3]) and L 2 ([0, 0.3]) errors of the approximate solution as a function 
of the synchronization time step, At. One can see a second-order dependence, i.e. the error scales 
like 0((At) 2 ). 

5.2. Simulation of a Diamond Network. An important point of particleclaw is that it tends 
to yield high quality numerical solutions (including shocks) for very few particles (see e.g. [TTlfT3] ). 
Here, we demonstrate this feature for a so-called diamond network. This network of seven edges 
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Figure 3. Bottleneck test case. Solution at time t — Figure 4. Bottleneck test case. The error as 
(top), t — 1.5 (middle), and t — 3 (bottom), computed a function of the synchronization time step At 
with particleclaw with h = 8 ■ 10 -2 and d = 2 • 10 -2 . exhibits second order convergence behavior. 

and six nodes is shown in Fig. [5] Vehicles enter in the node marked "A" driving towards "1"; 
then the traffic flow splits into a N and a E direction. At "2", drivers again decide to go SE or 
straight E. At "3" and "4" , two confluences happen, and finally vehicles exit the network at "B" . 
We choose all roads identical it™ = 1 an d = 1, and at all network nodes, an even split of traffic 
flow occurs. The lengths of all the roads is taken to be 1 (even if it does not seem so in the plot). 
The initial conditions are also identical on all the roads: u(x, 0) = 0.4 + 0.4cos(37r:r). We evolve 
the solution until the time ifi na i = 2. Figure [5] shows the computational results when applying 
particleclaw, with three different types of resolutions (in space and time). The top solution has 
an initial/desired distance of particles of h = 2 ■ 10 -2 , the distance of particles that are inserted 
near a shock is d = 5 • 10~ 3 , and the synchronization time step is At = 1 • 10~ 2 . In the middle 
solution, all of these numbers are multiplied by a factor of 5, and in the bottom solution by a 
factor of 20. That means that the top possesses 20 times the resolution than the bottom, and it is 
therefore about 400 times as costly to compute. One can observe that the higher resolution at the 
top yields sharper features. However, the presence and position of the shocks (and other features) 
is well captured by the resolution in the middle, and even at the bottom most features are visible 
quite clearly. 



6. Conclusions and Outlook 

The framework presented in this paper demonstrates that the characteristic particle method 
particleclaw can be generalized to solving first order traffic models on road networks in a robust, 
accurate, and exactly conservative fashion. The presented methodology allows to apply the basic 
one-dimensional solver on each edge. The coupling of traffic states at network nodes is achieved 
by a special synchronization step that ensures exact conservation properties by applying a proper 
transfer of area between edges. A fundamental ingredient in this synchronization is the invoking of 
generalized Riemann solvers that have been developed for traffic networks. In this paper, a specific 
focus lies on 1-to-l, l-to-2, and 2-to-l highway junctions. However, the presented methodology 
applies to more general types of networks, as long as the evolution on each edge is described by 
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Figure 5. Diamond network test case with various resolutions. The solutions have parameters h = 
2 • 10~ 2 s, d = 5 ■ 10 _3 s, and At = 1 ■ 10 _2 s, where s — 1 at the top, s — 5 in the middle, and s = 20 at 
the bottom. 

a scalar hyperbolic conservation law, and the problem description allows to formulate coupling 
conditions at the network nodes. 

As previously shown [TT], particleclaw itself can be made second-order accurate with respect 
to the spacing between particles. The numerical examples investigated here confirm that the 
presented method is also second order accurate with respect to the time step At between synchro- 
nization events. Moreover, it is demonstrated that the presented approach yields good quality 
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numerical solutions (with an accurate location of shocks) , even when using only few particles on 
each edge. 

In contrast to many traditional numerical methods, here the synchronization time step At is not 
limited by the spacing between particles (i.e. no CFL stability restriction) . It is only limited by the 
length of the edges, insofar as between two synchronization events, information is not allowed to 
travel beyond a single edge. In general, this allows for much larger values for At than traditional 
methods do. An exception is posed by networks that possess a few very short edges. These can be 
treated significantly more efficiently by a simple generalization of the presented approach, namely 
by using a different At for each network node, depending on the shortest edge connected to it. 
Thus, edges that neighbor short edges would need to synchronize more often, but the rest of the 
network would not. 

The accurate location of shocks with only few particles on each edge, and the possibility of 
invoking the coupling conditions at network nodes relatively rarely, allow for fast and memory- 
efficient implementations. This makes the approach attractive for the simulation of large road 
networks, as well as the optimization, design, and control of the traffic flow on such large networks. 
It is also conceivable that the presented methodology can be adapted to other kinds of nonlinear 
flows on networks, for instance continuum models of supply chains [Tl 120] . 
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